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We employ a recently developed model that allows the study of two-dimensional brittle crack 
propagation under fixed grip boundary conditions. The crack development highlights the importance 
of voids which appear ahead of the crack as observed in experiments on the nano-scale. The 
appearance of these voids is responsible for roughening the crack path on small scales, in agreement 
with theoretical expectations. With increasing speed of propagation one observes the branching 
instabilities that were reported in experiments. The simulations allow understanding the phenomena 
by analyzing the elastic stress field that accompanies the crack dynamics. 



I. INTRODUCTION 

In laboratory experiments one studies crack propaga- 
tion in elastic media by holding a stretched slab of ma- 
terial with a grip, and initiating the crack by making 
a small notch at one end of the sample pj. When the 
notch exceeds the Griffith's length (H, the crack propa- 
gates with increasing speed until there is a balance be- 
tween the release of elastic energy and the creation of 
surface energy [1, H| . It turned out that molecular dy- 
namics simulations failed systematically to mimic this set 
up, forcing simulators to pull continuously on the bound- 
aries to achieve a propagating crack With fixed grip 
boundary conditions crack tended to slow down and stop 
propagating. The inability to simulate a brittle crack 
by molecular dynamics gave rise to claims that it would 
be impossible to advance a brittle crack without contin- 
uous stretching. The fundamental reason for this long- 
standing problem was understood recently. In Ref. [(| it 
was shown that the solution of this conundrum lies in the 
range of the inter-particle potential. By changing the po- 
tential range one can go from a brittle to a ductile crack, 
and the latter is stopped by plastic dissipation. Only a 
brittle crack can support itself with a fixed grip. Brittlc- 
ness was shown to be guaranteed by choosing a potential 
range that is of the order of the inter-particle distance. 
In this paper we present results of molecular dynamics 
simulations of brittle cracks under a fixed grip that are 
based on this understanding. 

One of the interesting issues that can be studied by 
such simulations is the proposed existence of voids that 
open up before the crack, acting as pointers for further 
propagation. This issue had been quite controversial. 
While some experiments, and in particular those per- 
taining to ultra-slow crack propagation in glass in wet 
atmosphere, strongly indicated the appearance of voids 
ahead of the crack [7[ , other experiments failed to observe 
such voids, and claimed that they are irrelevant to the 
propagation of brittle cracks. In a theoretical study, Ref. 
Q constructed a model of void-dominated crack prop- 
agation and attempted to explain the observed rough- 
ness of brittle cracks by the randomness associated with 
the positioning of the voids ahead of the crack. The nu- 
merical simulations shown below appear to vindicate this 



approach almost verbatim. Voids do appear, and their 
appearance is random as stipulated in Ref. [Sj, and see 
Sect. Ofor more details. 

Another issue of interest is the instabilities of the prop- 
agating crack when its velocity increases. The simula- 
tions show very clearly branching instabilities with sub- 
sequent competition between the two branches that typ- 
ically result in the death of one branch in favor of the 
other. 

In Sect. [IQwc present the details of the preparation of 
the samples for the molecular dynamics simulations. In 
Sect. Mil we show that our cracks are brittle by demon- 
strating that they start running precisely when the Grif- 
fith's criterion met. In Sect. |lV]we discuss the appear- 
ance of voids ahead of the crack and the resulting rough- 
ening of the crack path. Finally, in Sect. [V]we present 
a brief discussion of the observed instabilities. The last 
section presents a short summary and indicates the road 
ahead. 



II. MOLECULAR DYNAMICS 
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FIG. 1: Color online: Examples of three different potentials 
with r co = 1.2, 1.8 and 2.4 generated using Eq. fl]). We used 
the potential with r co = 1.2 for the purpose of this paper. 



For the numerical experiments we employ a generic 
glass former in 2-dimcnsions in the form of a 50-50 binary 
mixture of 'small' and 'large' particles, chosen to avoid 
any crystallization. In fact the particles interact by inter- 
particle potentials as shown in Fig. [I] with the analytic 
form 
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Here r m i n /\ij is the length where the potential attain 
it's minimum, and r co /Xij is the cut-off length for which 
the potential vanishes. The coefficients a, b and cu are 
chosen such that the repulsive and attractive parts of 
the potential are continuous with two derivatives at the 
potential minimum and the potential goes to zero contin- 
uously at r co /Xij with two continuous derivatives as well. 
The interaction length-scale \j between any two parti- 
cles i and j is Xij — 1.0A, A,j = 1.18A and A,j = 1.4A for 
two 'small' particles, one 'large' and one 'small' particle 
and two 'large' particle respectively. The unit of length 
A is set to be the interaction length scale of two small 
particles, e is the unit of energy and ks = 1- 




FIG. 2: Color Online: a typical crack in a system of size 
1000Ax3000A which is self sustaining under fixed grip condi- 
tions with 7 = 2.5%. The initial cut was of length L = 750A. 
Note the roughening on small scales which is discussed in Sect. 

hyi 

We used a modified Berendsen thermostat which cou- 
ples a constant number of particles to the bath, regardless 
of the system size @. The preparation protocol starts 
with equilibrating the sample at high temperature and 
pressure, using periodic boundary conditions. The sys- 
tems are then cooled to temperature T = 10~ 2 while 
keeping the pressure high, in order to avoid the creation 
of any holes in the material. Then the pressure is reduced 
to zero, P = 0.0, such that the periodic boundary condi- 
tions could be removed; particles forming the right and 
left walls were frozen, but the upper and lower bound- 
aries were rendered free. The brittle crack experiment 
starts by loading the system uniaxially (with a constant 
velocity such that t> wa ii <C c s /10) until a desired stress 
is reached, and then the side walls are held fixed. A cut 
is then implemented by the cancelation of forces cross- 
ing an imaginary line of desired length which starts at 
the lower boundary. The evolution of the crack is simu- 
lated by molecular dynamics coupled to a heat bath at 
T = 0.0 and requires no further loading of the system. 
The results presented below were created using two dif- 
ferent geometries, one of width I000A and length 3000A, 



and the other being a square sample of I600A x 1600A. 
The density is the zero temperature and zero pressure 
density p = 0.745. 

In Fig. [5] we present a typical crack that results from 
this procedure, in which the loading 7 was 7 = 2.5% and 
an initial crack of length 750A. One observes the typical 
roughness that is discussed below in Sect. [IVJ 

III. THE GRIFFITH LENGTH 

Upon making the initial cut, the system always has a 
microscopic plastic response which however stops if the 
length of the cut is smaller than the Griffith length. To 
demonstrate that our running cracks are indeed brittle 
we test the Griffith's criterion in our simulation. 

The length of the initial cut that evolves spontaneously 
into a crack was determined by Griffith [2|, comparing 
the energy release from the stress field to the energy con- 
sumed by the creation of two new fresh surfaces by the 
advancing crack. With the stress exerted on the slab be- 
ing er, the critical length L c is determined by the bulk 
modulus E and the surface energy per unit length e via 
the relation 

In our simulations we are able to measure the critical 
length L c for a given loaded system. Since all the other 
material parameters appearing in Eq. [5] are measured 
independently, we can show that our model crack propa- 
gation agrees with the expected physics of brittle cracks. 

The stress field a is simply determined by our grip 
boundary conditions. The bulk modulus E is read from 
the linear part of a uniaxial straining experiment shown 
in Fig. |31 To estimate the surface energy in our system 
we computed the total energy before and after making 
the initial cut of length L. Taking the difference and 
dividing by L we get the estimate 

e « 15 . (3) 

At this point we employed 50 independent samples 
that were equilibrated at high temperature and quenched 
to the glassy state. The glassy samples were then 
loaded quasi-statically to strain values in the range 7 = 
1.6,1.8,2.0,2.2,2.4,2.6. In order to probe the Griffith's 
Length L c we tracked the behavior of our sample during 
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FIG. 3: Color Online: A typical stress-strain curve measured 
from a uniaxial straining experiment. We read the Bulk mod- 

° 1 1040 

ulus E ~ 725 from the slope and evaluate the stress for the 
theoretical estimates of the Griffith's Length. 
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FIG. 4: Color Online: Upper Panel: The linear relation- 
ship confirms the compliance of our model to Griffith's Law. 
The slope corresponds to the theoretical estimation emanat- 
ing from Eq. ((2| . Lower Panel: Plotted are the average values 
of the crack length C within 30 time units of the experiment. 
We read the critial Griffth's length L c at the size of L for 
which £ > 0. 



the first 30 time units of our simulation to see whether a 
crack has started to run as a result of the introduced cut. 
For each sample we repeated this procedure for various 
values of L. The results, (averaged over 50 samples) are 
presented in the bottom panel of Fig. |4] We show the 
length of the evolving crack, denoted as C, as a function 
the initial cut size L. Where the length C begins to run 
determines L c for the conditions at hand. 

In the upper panel of the same figure we plot a 2 vs. 
l/L c . The resulting linear plot is in agreement with Eq. 
([2|), with the slope being close to the theoretical expec- 
tation of 2Ee/ir. 



FIG. 5: Color Online: A typical crack tip with the voids that 
opened ahead of the crack 



IV. ROUGHENING RESULTING FROM VOIDS 
AHEAD OF THE CRACK 



The present numerical simulations allows the verifica- 
tion of the role of voids that form ahead of the crack, 
determining, at least on short length scales, the path cho- 
sen by the crack as it develops in the amorphous solid. 
We note that an amorphous solids is not an ideal elas- 
tic material, in which mathematically straight cracks are 
possible. The randomness of the material must show up 
in the geometry of the crack in this way or another. The 
voids ahead of the crack serve as pointers for the forth- 
coming propagation of the crack, and as is shown below, 
their random positions determines the roughening of the 
crack path. 

In Fig. [5] we show a typical crack tip with the voids 
that opened ahead. The physical reason why voids may 
prefer to open ahead of the crack and not on the crack 
tip was explained in Ref. . The argument is based on 
the fact that voids open due to plastic yield, and they do 
this where the pressure is maximal. Near the crack tip 
there is a process zone where the pressure is increasing 
going outward, until one hits the maximal pressure curve 
which connects with the outer elastic solution, see Fig. [5] 
In this figure we show the average pressure as measured 
in the numerical simulation in a fixed window tracking in 
time the crack tip. In particular one should pay attention 
to the lower right panel in Fig. [5] which shows that the 
maximal pressure lies ahead of the tip, at a distance £ 
from the tip. This figure should be compared with Fig. 2 
of Ref. which it directly vindicates. 

Obviously, the void that opens ahead of the tip can 
have a broadly distributed position around the forward 
direction (where the probability to form a void is max- 
imal). An actual measurement of this distribution for 
the present simulation is shown in Fig. [7] Shown are 
the actual positions of the voids as they appear ahead 
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FIG. 6: Color Online: Upper Paner: The pressure field in a 
fixed frame tracking the crack tip, averaged over time. The 
color code is brown for lowest pressure and blue the high- 
est. The two black lines are the horizontal and vertical cross 
section of this measured pressure field and are shown in the 
bottom left and right panels respectively. This last picture 
demonstrates the fact that the maximal pressure lies ahead 
of the crack tip where voids are being created. 
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FIG. 7: Color Online: Upper panel: the actual positions of the 
voids appearing ahead of the crack during the time evolution 
of the crack. These results are histogrammed in the lower two 
panels to provide the angular distribution of the void positions 
and the radial distribution in the forward direction. 



of the crack during the time evolution, and in the lower 
panels the probability distribution for the void to appear 
at an angle 9 with respect to the forward direction, and 
at a distance r along the forward direction. These data 
should be compared with Fig. 3 of Ref. §]■ 

Two points are worth stressing: first, the angular dis- 



tribution of void positions will be the source of roughen- 
ing of the crack - the probability to fall along the forward 
direction is not high enough compared to positive or neg- 
ative angles with the respect to the forward direction. 
Second, once a void appears on, say, a positive angle, 
the next void will have an even higher probability for a 
positive angle, meaning that the rough crack is expected 
to show a roughening exponent higher than 1/2. Indeed, 
such a persistent random walk is always expected to show 
exponents higher than 1/2, whereas anti-persistent ran- 
dom walks are characterized by a roughening exponent 
smaller than 1/2. 

The roughening exponent was measure here in the 
usual way, i.e. considering the crack as a graph y{x) 
and determining h(r) according to 

h(r) = (m&x{y(x}} x<x<x+r -mm{y(x} x<x<x+ r) x ■ (4) 

For a self-affine graph the scaling exponent C is defined 
via the scaling relation 



h(r) 



(5) 



We show this quantity in a log-log plot in Fig. [5] As ex- 
pected the function exhibits a persistent scaling exponent 
of C « 0.66 for scales r < 30, and then a cross over to a 
random graph without persistence or anti-persistence for 
higher scales. 



V. INSTABILITIES 

Crack dynamics are governed by the balance of energy 
at the crack tip. The influx of elastic energy through 
the crack-tip is used to create surface energy behind the 
advancing tip. With increasing crack velocity there is 
not sufficient surface to store the energy that is released. 
Therefore the system resorts to instabilities in the form 
of branching and oscillations in order to increase the 
amount of surface created per unit length on the axis of 
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FIG. 9: Color online: Upper panel: The crack length as a 
function of time; the velocity v is read from the slope. Lower 
panel: A typical branching event. One sees the two crack tips 
advancing simultaneously for a while, until the competition 
between them selects one or the other. Here the right crack 
wins, as is evidenced by the void ahead of it. 



propagation. In our simulations we find that when the ve- 



locity of the crack tip reaches about 30% of the Rayleigh 
speed, one begins to observe crack branching. We note 
that our system never develops two independently prop- 
agating cracks, but rather reduced the velocity of propa- 
gation through attempted branchings. Oscillations were 
not observed. A typical branching event is shown in Fig. 
[9J where we see the two crack tips as they still grow simul- 
taneously. The competition between them always results 
in the demise of one of the growing cracks in favor of the 
other, until the next branching event. 

VI. SUMMARY AND DISCUSSION 

In summary, we have demonstrated that molecular dy- 
namics can be usefully employed to study brittle crack 
propagation by selecting an appropriate interparticle po- 
tential. The resulting cracks are growing by nucleating 
voids ahead of the crack tip in much the same way that 
was anticipated by Ref. pj and theorized in Ref. 
This mechanism of growth is responsible for roughening 
of the crack path on small scales, but as long as side 
branching does not commence, the crack is glob ally flat 
on macroscopic scales, as expected theoretically |10| . 

In future work molecular dynamics simulations will be 
employed in three dimensions where the issue of micro- 
branching Q can be studied and compared with experi- 
mental results. 
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